# main screening function for genes
runScreening_oneGene <- function (gene) {
message ("Process gene " , gene, " (" , match (gene, genes), "/" , length (genes), ")" )
dat_gene <- dat_KO[dat_KO$ FeatureID == gene,]
# 1) main model estimation, based on model selection
estimate_ZImodels <- sum (dat_gene$ tpmAbd == 0 ) >= (0.1 * nrow (dat_gene))
if (estimate_ZImodels) { # estimate zero-inflated models
# these zero-inflated models sometimes run into an estimation error. In this case, ignore the respective model.
# If both models run into an estimation error, simply estimate the non-zero-inflated models instead.
model_list <- list ("ZI Gaussian" = estimate_model (dat_model = dat_gene,
model_type = "ZI Gaussian" ,
model_formula = as.formula ("tpmAbd ~ intervention + phase + gender + age_centered50 + bmi_centered25" )),
"ZI log-Gaussian" = estimate_model (dat_model = dat_gene,
model_type = "ZI Gaussian" ,
model_formula = as.formula ("tpmAbd_log ~ intervention + phase + gender + age_centered50 + bmi_centered25" )))
fallBack_toNonZImodels <- ifelse (is.null (model_list[[1 ]]) & is.null (model_list[[2 ]]), TRUE , FALSE )
if (! fallBack_toNonZImodels) {
if (is.null (model_list[[1 ]]) | is.null (model_list[[2 ]])) {
null_index <- unname (which (sapply (model_list, function (x) { is.null (x) })))
model_list <- model_list[- null_index]
}
bestModel_index <- unname (which.min (sapply (model_list, function (x) { mean (abs (summary (x)$ residuals)) })))
model <- model_list[[bestModel_index]]
model_label <- names (model_list)[bestModel_index]
tab <- summary (model)$ tTable
colnames (tab) <- c ("Estimate" , "Std. Error" , "DF" , "t-value" , "Pr(>|t|)" ) # make names consistent to gam tab names
}
} else {
fallBack_toNonZImodels <- FALSE # simply to have this object defined also if ZI models are not estimated
}
if (! estimate_ZImodels || fallBack_toNonZImodels) { # estimate models without zero-inflation
model_list <- list ("Gaussian" = estimate_model (dat_model = dat_gene,
model_type = "Gaussian" ,
model_formula = as.formula ("tpmAbd ~ intervention + phase + gender + age_centered50 + bmi_centered25 + s(participant_id, bs = 're')" )),
"log-Gaussian" = estimate_model (dat_model = dat_gene,
model_type = "Gaussian" ,
model_formula = as.formula ("tpmAbd_log ~ intervention + phase + gender + age_centered50 + bmi_centered25 + s(participant_id, bs = 're')" )))
bestModel_index <- unname (which.max (sapply (model_list, function (x) { summary (x)$ dev.expl })))
model <- model_list[[bestModel_index]]
model_label <- names (model_list)[bestModel_index]
tab <- summary (model)$ p.table
}
# 2) interaction model estimation, only for the final selected model
# model estimation with the lower category of the interaction variable as the reference category
model_intList <- estimate_interactionModels (dat_model = dat_gene,
selected_model = model_label)
# model estimation with the higher category of the interaction variable as the reference category
dat_gene_changedRef <- dat_gene
dat_gene_changedRef$ gender <- relevel (dat_gene_changedRef$ gender, ref = "W" )
dat_gene_changedRef$ age_cat <- relevel (dat_gene_changedRef$ age_cat, ref = "50-69" )
dat_gene_changedRef$ bmi_cat <- relevel (dat_gene_changedRef$ bmi_cat, ref = "[25, 31]" )
dat_gene_changedRef$ dailyFiber_cat <- relevel (dat_gene_changedRef$ dailyFiber_cat, ref = ">= 30g" )
dat_gene_changedRef$ baselineDiv_cat <- relevel (dat_gene_changedRef$ baselineDiv_cat, ref = "[3.85, 4.4)" )
model_intList_changedRef <- estimate_interactionModels (dat_model = dat_gene_changedRef,
selected_model = model_label)
# create joint table
tab_intEffects <- data.frame (effect_type = c ("fresh intervention" , "pasteurized intervention" ),
effect_ageLow = c (model_intList$ tab_ageInt["interventionFresh" , "Estimate" ],
model_intList$ tab_ageInt["interventionPasteurized" , "Estimate" ]),
effect_ageHigh = c (model_intList_changedRef$ tab_ageInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_ageInt["interventionPasteurized" , "Estimate" ]),
se_ageLow = c (model_intList$ tab_ageInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_ageInt["interventionPasteurized" , "Std. Error" ]),
se_ageHigh = c (model_intList_changedRef$ tab_ageInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_ageInt["interventionPasteurized" , "Std. Error" ]),
pvalue_ageLow = c (model_intList$ tab_ageInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_ageInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_ageHigh = c (model_intList_changedRef$ tab_ageInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_ageInt["interventionPasteurized" , "Pr(>|t|)" ]),
effect_sexM = c (model_intList$ tab_sexInt["interventionFresh" , "Estimate" ],
model_intList$ tab_sexInt["interventionPasteurized" , "Estimate" ]),
effect_sexW = c (model_intList_changedRef$ tab_sexInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_sexInt["interventionPasteurized" , "Estimate" ]),
se_sexM = c (model_intList$ tab_sexInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_sexInt["interventionPasteurized" , "Std. Error" ]),
se_sexW = c (model_intList_changedRef$ tab_sexInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_sexInt["interventionPasteurized" , "Std. Error" ]),
pvalue_sexM = c (model_intList$ tab_sexInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_sexInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_sexW = c (model_intList_changedRef$ tab_sexInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_sexInt["interventionPasteurized" , "Pr(>|t|)" ]),
effect_bmiLow = c (model_intList$ tab_bmiInt["interventionFresh" , "Estimate" ],
model_intList$ tab_bmiInt["interventionPasteurized" , "Estimate" ]),
effect_bmiHigh = c (model_intList_changedRef$ tab_bmiInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_bmiInt["interventionPasteurized" , "Estimate" ]),
se_bmiLow = c (model_intList$ tab_bmiInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_bmiInt["interventionPasteurized" , "Std. Error" ]),
se_bmiHigh = c (model_intList_changedRef$ tab_bmiInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_bmiInt["interventionPasteurized" , "Std. Error" ]),
pvalue_bmiLow = c (model_intList$ tab_bmiInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_bmiInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_bmiHigh = c (model_intList_changedRef$ tab_bmiInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_bmiInt["interventionPasteurized" , "Pr(>|t|)" ]),
effect_fibLow = c (model_intList$ tab_fibInt["interventionFresh" , "Estimate" ],
model_intList$ tab_fibInt["interventionPasteurized" , "Estimate" ]),
effect_fibHigh = c (model_intList_changedRef$ tab_fibInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_fibInt["interventionPasteurized" , "Estimate" ]),
se_fibLow = c (model_intList$ tab_fibInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_fibInt["interventionPasteurized" , "Std. Error" ]),
se_fibHigh = c (model_intList_changedRef$ tab_fibInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_fibInt["interventionPasteurized" , "Std. Error" ]),
pvalue_fibLow = c (model_intList$ tab_fibInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_fibInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_fibHigh = c (model_intList_changedRef$ tab_fibInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_fibInt["interventionPasteurized" , "Pr(>|t|)" ]),
effect_divLow = c (model_intList$ tab_divInt["interventionFresh" , "Estimate" ],
model_intList$ tab_divInt["interventionPasteurized" , "Estimate" ]),
effect_divHigh = c (model_intList_changedRef$ tab_divInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_divInt["interventionPasteurized" , "Estimate" ]),
se_divLow = c (model_intList$ tab_divInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_divInt["interventionPasteurized" , "Std. Error" ]),
se_divHigh = c (model_intList_changedRef$ tab_divInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_divInt["interventionPasteurized" , "Std. Error" ]),
pvalue_divLow = c (model_intList$ tab_divInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_divInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_divHigh = c (model_intList_changedRef$ tab_divInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_divInt["interventionPasteurized" , "Pr(>|t|)" ]))
# 3) gather results
data.frame (gene = gene,
relFreq_zeroTPM = sum (dat_gene$ tpmAbd == 0 ) / nrow (dat_gene),
relFreq_min10TPM = sum (dat_gene$ tpmAbd >= 10 ) / nrow (dat_gene),
min_tpmAbd = min (dat_gene$ tpmAbd),
mean_tpmAbd = mean (dat_gene$ tpmAbd),
median_tpmAbd = median (dat_gene$ tpmAbd),
max_tpmAbd = max (dat_gene$ tpmAbd),
baseline_medianTPMAbd = median (dat_gene$ tpmAbd[dat_gene$ inter_treat %in% c ("Baseline.Fresh" ,"Baseline.Pasteurized" )]),
baseline_prevalence = unname (prop.table (table (dat_gene$ tpmAbd[dat_gene$ inter_treat %in% c ("Baseline.Fresh" ,"Baseline.Pasteurized" )] > 10 ))["TRUE" ]),
postFresh_medianTPMAbd = median (dat_gene$ tpmAbd[dat_gene$ inter_treat == "After.Fresh" ]),
postFresh_prevalence = unname (prop.table (table (dat_gene$ tpmAbd[dat_gene$ inter_treat == "After.Fresh" ] > 10 ))["TRUE" ]),
postPast_medianTPMAbd = median (dat_gene$ tpmAbd[dat_gene$ inter_treat == "After.Pasteurized" ]),
postPast_prevalence = unname (prop.table (table (dat_gene$ tpmAbd[dat_gene$ inter_treat == "After.Pasteurized" ] > 10 ))["TRUE" ]),
model = model_label,
model_deviance = ifelse (! grepl ("ZI" , model_label), summary (model)$ dev.expl, NA ),
model_meanAbsResid = ifelse (grepl ("ZI" , model_label), mean (abs (summary (model)$ residuals)), NA ),
freshInt_effect = tab["interventionFresh" , "Estimate" ],
freshInt_se = tab["interventionFresh" , "Std. Error" ],
freshInt_pvalue = tab["interventionFresh" , "Pr(>|t|)" ],
pastInt_effect = tab["interventionPasteurized" , "Estimate" ],
pastInt_se = tab["interventionPasteurized" , "Std. Error" ],
pastInt_pvalue = tab["interventionPasteurized" , "Pr(>|t|)" ],
phaseEffect = tab["phasephase 2" , "Estimate" ],
genderEffect = tab["genderW" , "Estimate" ],
ageEffect = tab["age_centered50" , "Estimate" ],
bmiEffect = tab["bmi_centered25" , "Estimate" ],
freshInt_ageLow_effect = tab_intEffects$ effect_ageLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageHigh_effect = tab_intEffects$ effect_ageHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageLow_se = tab_intEffects$ se_ageLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageHigh_se = tab_intEffects$ se_ageHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageLow_pvalue = tab_intEffects$ pvalue_ageLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageHigh_pvalue = tab_intEffects$ pvalue_ageHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexM_effect = tab_intEffects$ effect_sexM[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexW_effect = tab_intEffects$ effect_sexW[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexM_se = tab_intEffects$ se_sexM[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexW_se = tab_intEffects$ se_sexW[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexM_pvalue = tab_intEffects$ pvalue_sexM[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexW_pvalue = tab_intEffects$ pvalue_sexW[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiLow_effect = tab_intEffects$ effect_bmiLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiHigh_effect = tab_intEffects$ effect_bmiHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiLow_se = tab_intEffects$ se_bmiLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiHigh_se = tab_intEffects$ se_bmiHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiLow_pvalue = tab_intEffects$ pvalue_bmiLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiHigh_pvalue = tab_intEffects$ pvalue_bmiHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibLow_effect = tab_intEffects$ effect_fibLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibHigh_effect = tab_intEffects$ effect_fibHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibLow_se = tab_intEffects$ se_fibLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibHigh_se = tab_intEffects$ se_fibHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibLow_pvalue = tab_intEffects$ pvalue_fibLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibHigh_pvalue = tab_intEffects$ pvalue_fibHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divLow_effect = tab_intEffects$ effect_divLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divHigh_effect = tab_intEffects$ effect_divHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divLow_se = tab_intEffects$ se_divLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divHigh_se = tab_intEffects$ se_divHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divLow_pvalue = tab_intEffects$ pvalue_divLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divHigh_pvalue = tab_intEffects$ pvalue_divHigh[tab_intEffects$ effect_type == "fresh intervention" ],
pastInt_ageLow_effect = tab_intEffects$ effect_ageLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageHigh_effect = tab_intEffects$ effect_ageHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageLow_se = tab_intEffects$ se_ageLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageHigh_se = tab_intEffects$ se_ageHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageLow_pvalue = tab_intEffects$ pvalue_ageLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageHigh_pvalue = tab_intEffects$ pvalue_ageHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexM_effect = tab_intEffects$ effect_sexM[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexW_effect = tab_intEffects$ effect_sexW[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexM_se = tab_intEffects$ se_sexM[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexW_se = tab_intEffects$ se_sexW[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexM_pvalue = tab_intEffects$ pvalue_sexM[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexW_pvalue = tab_intEffects$ pvalue_sexW[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiLow_effect = tab_intEffects$ effect_bmiLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiHigh_effect = tab_intEffects$ effect_bmiHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiLow_se = tab_intEffects$ se_bmiLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiHigh_se = tab_intEffects$ se_bmiHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiLow_pvalue = tab_intEffects$ pvalue_bmiLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiHigh_pvalue = tab_intEffects$ pvalue_bmiHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibLow_effect = tab_intEffects$ effect_fibLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibHigh_effect = tab_intEffects$ effect_fibHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibLow_se = tab_intEffects$ se_fibLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibHigh_se = tab_intEffects$ se_fibHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibLow_pvalue = tab_intEffects$ pvalue_fibLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibHigh_pvalue = tab_intEffects$ pvalue_fibHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divLow_effect = tab_intEffects$ effect_divLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divHigh_effect = tab_intEffects$ effect_divHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divLow_se = tab_intEffects$ se_divLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divHigh_se = tab_intEffects$ se_divHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divLow_pvalue = tab_intEffects$ pvalue_divLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divHigh_pvalue = tab_intEffects$ pvalue_divHigh[tab_intEffects$ effect_type == "pasteurized intervention" ])
}
# main screening function for pathways
runScreening_onePathway <- function (pathway) {
message ("Process pathway " , pathway, " (" , match (pathway, pathways), "/" , length (pathways), ")" )
dat_pathway <- dat_PW[dat_PW$ taxonomy_level3 == pathway,]
# 1) main model estimation, based on model selection
estimate_ZImodels <- sum (dat_pathway$ tpmAbd == 0 ) >= (0.1 * nrow (dat_pathway))
if (estimate_ZImodels) { # estimate zero-inflated models
# these zero-inflated models sometimes run into an estimation error. In this case, ignore the respective model.
# If both models run into an estimation error, simply estimate the non-zero-inflated models instead.
model_list <- list ("ZI Gaussian" = estimate_model (dat_model = dat_pathway,
model_type = "ZI Gaussian" ,
model_formula = as.formula ("tpmAbd ~ intervention + phase + gender + age_centered50 + bmi_centered25" )),
"ZI log-Gaussian" = estimate_model (dat_model = dat_pathway,
model_type = "ZI Gaussian" ,
model_formula = as.formula ("tpmAbd_log ~ intervention + phase + gender + age_centered50 + bmi_centered25" )))
fallBack_toNonZImodels <- ifelse (is.null (model_list[[1 ]]) & is.null (model_list[[2 ]]), TRUE , FALSE )
if (! fallBack_toNonZImodels) {
if (is.null (model_list[[1 ]]) | is.null (model_list[[2 ]])) {
null_index <- unname (which (sapply (model_list, function (x) { is.null (x) })))
model_list <- model_list[- null_index]
}
bestModel_index <- unname (which.min (sapply (model_list, function (x) { mean (abs (summary (x)$ residuals)) })))
model <- model_list[[bestModel_index]]
model_label <- names (model_list)[bestModel_index]
tab <- summary (model)$ tTable
colnames (tab) <- c ("Estimate" , "Std. Error" , "DF" , "t-value" , "Pr(>|t|)" ) # make names consistent to gam tab names
}
} else {
fallBack_toNonZImodels <- FALSE # simply to have this object defined also if ZI models are not estimated
}
if (! estimate_ZImodels || fallBack_toNonZImodels) { # estimate models without zero-inflation
model_list <- list ("Gaussian" = estimate_model (dat_model = dat_pathway,
model_type = "Gaussian" ,
model_formula = as.formula ("tpmAbd ~ intervention + phase + gender + age_centered50 + bmi_centered25 + s(participant_id, bs = 're')" )),
"log-Gaussian" = estimate_model (dat_model = dat_pathway,
model_type = "Gaussian" ,
model_formula = as.formula ("tpmAbd_log ~ intervention + phase + gender + age_centered50 + bmi_centered25 + s(participant_id, bs = 're')" )))
bestModel_index <- unname (which.max (sapply (model_list, function (x) { summary (x)$ dev.expl })))
model <- model_list[[bestModel_index]]
model_label <- names (model_list)[bestModel_index]
tab <- summary (model)$ p.table
}
# 2) interaction model estimation, only for the final selected model
# model estimation with the lower category of the interaction variable as the reference category
model_intList <- estimate_interactionModels (dat_model = dat_pathway,
selected_model = model_label)
# model estimation with the higher category of the interaction variable as the reference category
dat_pathway_changedRef <- dat_pathway
dat_pathway_changedRef$ gender <- relevel (dat_pathway_changedRef$ gender, ref = "W" )
dat_pathway_changedRef$ age_cat <- relevel (dat_pathway_changedRef$ age_cat, ref = "50-69" )
dat_pathway_changedRef$ bmi_cat <- relevel (dat_pathway_changedRef$ bmi_cat, ref = "[25, 31]" )
dat_pathway_changedRef$ dailyFiber_cat <- relevel (dat_pathway_changedRef$ dailyFiber_cat, ref = ">= 30g" )
dat_pathway_changedRef$ baselineDiv_cat <- relevel (dat_pathway_changedRef$ baselineDiv_cat, ref = "[3.85, 4.4)" )
model_intList_changedRef <- estimate_interactionModels (dat_model = dat_pathway_changedRef,
selected_model = model_label)
# create joint table
tab_intEffects <- data.frame (effect_type = c ("fresh intervention" , "pasteurized intervention" ),
effect_ageLow = c (model_intList$ tab_ageInt["interventionFresh" , "Estimate" ],
model_intList$ tab_ageInt["interventionPasteurized" , "Estimate" ]),
effect_ageHigh = c (model_intList_changedRef$ tab_ageInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_ageInt["interventionPasteurized" , "Estimate" ]),
se_ageLow = c (model_intList$ tab_ageInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_ageInt["interventionPasteurized" , "Std. Error" ]),
se_ageHigh = c (model_intList_changedRef$ tab_ageInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_ageInt["interventionPasteurized" , "Std. Error" ]),
pvalue_ageLow = c (model_intList$ tab_ageInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_ageInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_ageHigh = c (model_intList_changedRef$ tab_ageInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_ageInt["interventionPasteurized" , "Pr(>|t|)" ]),
effect_sexM = c (model_intList$ tab_sexInt["interventionFresh" , "Estimate" ],
model_intList$ tab_sexInt["interventionPasteurized" , "Estimate" ]),
effect_sexW = c (model_intList_changedRef$ tab_sexInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_sexInt["interventionPasteurized" , "Estimate" ]),
se_sexM = c (model_intList$ tab_sexInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_sexInt["interventionPasteurized" , "Std. Error" ]),
se_sexW = c (model_intList_changedRef$ tab_sexInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_sexInt["interventionPasteurized" , "Std. Error" ]),
pvalue_sexM = c (model_intList$ tab_sexInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_sexInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_sexW = c (model_intList_changedRef$ tab_sexInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_sexInt["interventionPasteurized" , "Pr(>|t|)" ]),
effect_bmiLow = c (model_intList$ tab_bmiInt["interventionFresh" , "Estimate" ],
model_intList$ tab_bmiInt["interventionPasteurized" , "Estimate" ]),
effect_bmiHigh = c (model_intList_changedRef$ tab_bmiInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_bmiInt["interventionPasteurized" , "Estimate" ]),
se_bmiLow = c (model_intList$ tab_bmiInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_bmiInt["interventionPasteurized" , "Std. Error" ]),
se_bmiHigh = c (model_intList_changedRef$ tab_bmiInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_bmiInt["interventionPasteurized" , "Std. Error" ]),
pvalue_bmiLow = c (model_intList$ tab_bmiInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_bmiInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_bmiHigh = c (model_intList_changedRef$ tab_bmiInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_bmiInt["interventionPasteurized" , "Pr(>|t|)" ]),
effect_fibLow = c (model_intList$ tab_fibInt["interventionFresh" , "Estimate" ],
model_intList$ tab_fibInt["interventionPasteurized" , "Estimate" ]),
effect_fibHigh = c (model_intList_changedRef$ tab_fibInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_fibInt["interventionPasteurized" , "Estimate" ]),
se_fibLow = c (model_intList$ tab_fibInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_fibInt["interventionPasteurized" , "Std. Error" ]),
se_fibHigh = c (model_intList_changedRef$ tab_fibInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_fibInt["interventionPasteurized" , "Std. Error" ]),
pvalue_fibLow = c (model_intList$ tab_fibInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_fibInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_fibHigh = c (model_intList_changedRef$ tab_fibInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_fibInt["interventionPasteurized" , "Pr(>|t|)" ]),
effect_divLow = c (model_intList$ tab_divInt["interventionFresh" , "Estimate" ],
model_intList$ tab_divInt["interventionPasteurized" , "Estimate" ]),
effect_divHigh = c (model_intList_changedRef$ tab_divInt["interventionFresh" , "Estimate" ],
model_intList_changedRef$ tab_divInt["interventionPasteurized" , "Estimate" ]),
se_divLow = c (model_intList$ tab_divInt["interventionFresh" , "Std. Error" ],
model_intList$ tab_divInt["interventionPasteurized" , "Std. Error" ]),
se_divHigh = c (model_intList_changedRef$ tab_divInt["interventionFresh" , "Std. Error" ],
model_intList_changedRef$ tab_divInt["interventionPasteurized" , "Std. Error" ]),
pvalue_divLow = c (model_intList$ tab_divInt["interventionFresh" , "Pr(>|t|)" ],
model_intList$ tab_divInt["interventionPasteurized" , "Pr(>|t|)" ]),
pvalue_divHigh = c (model_intList_changedRef$ tab_divInt["interventionFresh" , "Pr(>|t|)" ],
model_intList_changedRef$ tab_divInt["interventionPasteurized" , "Pr(>|t|)" ]))
# 3) gather results
data.frame (pathway = pathway,
relFreq_zeroTPM = sum (dat_pathway$ tpmAbd == 0 ) / nrow (dat_pathway),
relFreq_min10TPM = sum (dat_pathway$ tpmAbd >= 10 ) / nrow (dat_pathway),
min_tpmAbd = min (dat_pathway$ tpmAbd),
mean_tpmAbd = mean (dat_pathway$ tpmAbd),
median_tpmAbd = median (dat_pathway$ tpmAbd),
max_tpmAbd = max (dat_pathway$ tpmAbd),
baseline_medianTPMAbd = median (dat_pathway$ tpmAbd[dat_pathway$ inter_treat %in% c ("Baseline.Fresh" ,"Baseline.Pasteurized" )]),
baseline_prevalence = unname (prop.table (table (dat_pathway$ tpmAbd[dat_pathway$ inter_treat %in% c ("Baseline.Fresh" ,"Baseline.Pasteurized" )] > 10 ))["TRUE" ]),
postFresh_medianTPMAbd = median (dat_pathway$ tpmAbd[dat_pathway$ inter_treat == "After.Fresh" ]),
postFresh_prevalence = unname (prop.table (table (dat_pathway$ tpmAbd[dat_pathway$ inter_treat == "After.Fresh" ] > 10 ))["TRUE" ]),
postPast_medianTPMAbd = median (dat_pathway$ tpmAbd[dat_pathway$ inter_treat == "After.Pasteurized" ]),
postPast_prevalence = unname (prop.table (table (dat_pathway$ tpmAbd[dat_pathway$ inter_treat == "After.Pasteurized" ] > 10 ))["TRUE" ]),
model = model_label,
model_deviance = ifelse (! grepl ("ZI" , model_label), summary (model)$ dev.expl, NA ),
model_meanAbsResid = ifelse (grepl ("ZI" , model_label), mean (abs (summary (model)$ residuals)), NA ),
freshInt_effect = tab["interventionFresh" , "Estimate" ],
freshInt_se = tab["interventionFresh" , "Std. Error" ],
freshInt_pvalue = tab["interventionFresh" , "Pr(>|t|)" ],
pastInt_effect = tab["interventionPasteurized" , "Estimate" ],
pastInt_se = tab["interventionPasteurized" , "Std. Error" ],
pastInt_pvalue = tab["interventionPasteurized" , "Pr(>|t|)" ],
phaseEffect = tab["phasephase 2" , "Estimate" ],
genderEffect = tab["genderW" , "Estimate" ],
ageEffect = tab["age_centered50" , "Estimate" ],
bmiEffect = tab["bmi_centered25" , "Estimate" ],
freshInt_ageLow_effect = tab_intEffects$ effect_ageLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageHigh_effect = tab_intEffects$ effect_ageHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageLow_se = tab_intEffects$ se_ageLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageHigh_se = tab_intEffects$ se_ageHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageLow_pvalue = tab_intEffects$ pvalue_ageLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_ageHigh_pvalue = tab_intEffects$ pvalue_ageHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexM_effect = tab_intEffects$ effect_sexM[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexW_effect = tab_intEffects$ effect_sexW[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexM_se = tab_intEffects$ se_sexM[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexW_se = tab_intEffects$ se_sexW[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexM_pvalue = tab_intEffects$ pvalue_sexM[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_sexW_pvalue = tab_intEffects$ pvalue_sexW[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiLow_effect = tab_intEffects$ effect_bmiLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiHigh_effect = tab_intEffects$ effect_bmiHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiLow_se = tab_intEffects$ se_bmiLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiHigh_se = tab_intEffects$ se_bmiHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiLow_pvalue = tab_intEffects$ pvalue_bmiLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_bmiHigh_pvalue = tab_intEffects$ pvalue_bmiHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibLow_effect = tab_intEffects$ effect_fibLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibHigh_effect = tab_intEffects$ effect_fibHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibLow_se = tab_intEffects$ se_fibLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibHigh_se = tab_intEffects$ se_fibHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibLow_pvalue = tab_intEffects$ pvalue_fibLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_fibHigh_pvalue = tab_intEffects$ pvalue_fibHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divLow_effect = tab_intEffects$ effect_divLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divHigh_effect = tab_intEffects$ effect_divHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divLow_se = tab_intEffects$ se_divLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divHigh_se = tab_intEffects$ se_divHigh[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divLow_pvalue = tab_intEffects$ pvalue_divLow[tab_intEffects$ effect_type == "fresh intervention" ],
freshInt_divHigh_pvalue = tab_intEffects$ pvalue_divHigh[tab_intEffects$ effect_type == "fresh intervention" ],
pastInt_ageLow_effect = tab_intEffects$ effect_ageLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageHigh_effect = tab_intEffects$ effect_ageHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageLow_se = tab_intEffects$ se_ageLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageHigh_se = tab_intEffects$ se_ageHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageLow_pvalue = tab_intEffects$ pvalue_ageLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_ageHigh_pvalue = tab_intEffects$ pvalue_ageHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexM_effect = tab_intEffects$ effect_sexM[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexW_effect = tab_intEffects$ effect_sexW[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexM_se = tab_intEffects$ se_sexM[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexW_se = tab_intEffects$ se_sexW[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexM_pvalue = tab_intEffects$ pvalue_sexM[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_sexW_pvalue = tab_intEffects$ pvalue_sexW[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiLow_effect = tab_intEffects$ effect_bmiLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiHigh_effect = tab_intEffects$ effect_bmiHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiLow_se = tab_intEffects$ se_bmiLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiHigh_se = tab_intEffects$ se_bmiHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiLow_pvalue = tab_intEffects$ pvalue_bmiLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_bmiHigh_pvalue = tab_intEffects$ pvalue_bmiHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibLow_effect = tab_intEffects$ effect_fibLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibHigh_effect = tab_intEffects$ effect_fibHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibLow_se = tab_intEffects$ se_fibLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibHigh_se = tab_intEffects$ se_fibHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibLow_pvalue = tab_intEffects$ pvalue_fibLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_fibHigh_pvalue = tab_intEffects$ pvalue_fibHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divLow_effect = tab_intEffects$ effect_divLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divHigh_effect = tab_intEffects$ effect_divHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divLow_se = tab_intEffects$ se_divLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divHigh_se = tab_intEffects$ se_divHigh[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divLow_pvalue = tab_intEffects$ pvalue_divLow[tab_intEffects$ effect_type == "pasteurized intervention" ],
pastInt_divHigh_pvalue = tab_intEffects$ pvalue_divHigh[tab_intEffects$ effect_type == "pasteurized intervention" ])
}